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ABSTRACT 

We investigate the structure of cold dark matter halos using advanced models of spherical collapse 
and accretion in an expanding Universe. These base on solving time-dependent equations for the 
moments of the phase-space distribution function in the fluid approximation; our approach includes 
non-radial random motions, and most importantly, an advanced treatment of both dynamical relax- 
ation effects that takes place in the infalling matter: phase-mixing associated to shell crossing, and 
collective collisions related to physical dumpiness. We find self-similar solutions for the spherically- 
averaged profiles of mass density p{r), pseudo phase-space density Q{r) and anisotropy parameter 
f3{r). These profiles agree with the outcomes of state-of-the-art iV— body simulations in the radial 
range currently probed by the latter; at smaller radii, we provide specific predictions. In the perspec- 
tive provided by our self-similar solutions we link the halo structure to its two-stage growth history, 
and propose the following picture. During the early fast collapse of the inner region dominated by a 
few merging clumps, efficient dynamical relaxation plays a key role in producing a closely universal 
mass density and pseudo phase-space density profiles; in particular, these are found to depend only 
weakly on the detailed shape of the initial perturbation and the related collapse times. The subsequent 
inside-out growth of the outer regions feeds on the slow accretion of many small clumps and diffuse 
matter; thus the outskirts are only mildly affected by dynamical relaxation but are more sensitive to 
asymmetries and cosmological variance. 

Subject headings: dark matter — galaxies: halos — method: analytical. 



1. INTRODUCTION 

The cosmogonical paradigm envisages galaxies and 
galaxy systems to form and shine when baryons accrete, 
settle and condense under the gravitational pull provided 
by the dark matter (DM) , which is cold and collisionless 
at the binary level. Understanding the formation history 
and the detailed structure of the DM gravitational wells 
is needed to provide a firm foundation for an accurate 
theory of galaxy formation and evolution. 

The story starts with DM primordial density pertur- 
bations of the cosmic density field; these at first grow 
by gravitational instability, and as the local gravity pre- 
vails are enforced to collapse and virialize into equilib- 
rium 'halos'. The resulting halo growth is hierarchical in 
mass and sequential in time, with small clumps forming 
first and then stochastically merging together into larger 
and more massive objects. 

However, describing and understanding the details of 
the halo formation process has proven to be a com- 
plex and tricky task, that still challenges the intense ef- 
forts spent by the astrophysical community over the last 
thirty years. Major empirical improvements occurred 
with the advent of the technology allowing to run inten- 
sive iV-body simulations of many million to a few billion 
particles on powerful supercomputers (for a review, see 
Springel et al. 2005). A number of hot issues are recalled 
next. 

First, DM halos form in two stages (Zhao et al. 2003; 
Diemand et al. 2007; Stadel et al. 2009; Fakhouri, Ma 
& Boylan-Kolchin 2010; Genel et al. 2010; Wang et al. 
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2011): an early fast collapse of the halo bulk including 
a few major merger events, which reshuffie the gravita- 
tional potential and cause the DM to undergo (incom- 
plete) dynamical relaxation; a late slow growth of the 
halo outskirts in the form of many minor mergers and 
diffuse accretion, which little affect the inner potential 
well, but contribute most (up to 80%) of the final mass. 

Second, DM halos in equilibrium show an approxi- 
mately universal spherically-averaged mass distribution. 
Originally, this was described with the simple formula 
proposed by NEW (Navarro, Frenk & White 1997), 
where the logarithmic density slope 7 = — dlogp/dlogr 
is given by 7 = (1 -I- 3 f)/(l -I- f) in terms of the radius 
f = r/r_2 normalized to the position r„2 where 7 = 2. 
This implies that the density profile p{r) oc (1 -|-f)~^ 



features an inner powerlaw behavior 



wards to 



steepens out- 



in the halo middle, and then goes into 



an asymptotic shape r~ in the outskirts. Remarkably, 
the NFW formula was found to be approximately scale- 
invariant, i.e., to fit accurately the simulation outcomes 
for halos of different mass scales. On the other hand, a 
weak scale dependence is introduced by the concentra- 
tion parameter c = i?20o/''-2, that is the ratio between 
r_2 and the radius i?2oo where the overdensity relative 
to the background amounts to 200; despite the name, 
it actually constitutes a measure of the halo outer ex- 
tension. Simulations show that c « 3.5 applies at the 
end of the fast collapse (with minor mass and redshift 
dependencies, see Prada et al. 2011), and increases as 
c oc H(z) oc (1 -|- z)~^ afterwards. At the present time 
z ~ 0, a relation c oc M^'^-^^ applies, with a scatter 
around 0.2 dex due to variance in the growth histories 
(see Bullock et al. 2001); typically, a galaxy halo with 
current mass M « a few 10^^ Mq collapsed at 2 « 2 
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features a concentration value c ~ 10, while the halo of 
a galaxy cluster with M k, Mq collapsed at ^ ~ 0.5 
features a value c ~ 5. 

In addition, recent simulations (Navarro et al. 2004, 
2010) have shown evidence of small but systematic de- 
viations from the simple NFW expression; they rather 
favor a Scrsic-Einasto (Scrsic 1963, Einasto 1965; see 
also Prugniel & Simien 1997; Merritt ct al. 2006; Lapi 
& Cavaliere 2011) formula 7 = r + (2 — r) f^', depend- 
ing on the two parameters r and rj that describe the 
inner behavior and the curvature of the density profile 
p{r) oc r""^ p-i^~'^)r'^ /''I before a final exponential cutoff. 
Current simulations set the upper limit r < 0.9 on the 
inner asymptotic powerlaw but still lack enough resolu- 
tion to pinpoint the true value. On the other hand, if a 
pure Einasto (1965) model with r = is adopted, dif- 
ferent values of 77 « 0.1 — 0.2 are required to provide 
precise fits to different simulated halos. This indicates 
that scale-invariancc in the mass distribution is actually 
broken, and/or that a halo's development is appreciably 
affected by variance related to detailed merging histories 
or environmental conditions. 

Third, DM halos show a definite spherically- aver aged 
profile of anisotropy. This is described on using the stan- 
dard Binncy (1978) parameter (3{r) = 1 — a0{r)/a^{r) 
in terms of the radial and tangential velocity disper- 
sions cr^ g{r). Halos are found to be quasi-isotropic at 
the center with /3 w within the resolution limits, and 
to become radially anisotropic outwards. Specifically, 
13 S3 0.25 applies at r « r_2 and /3 w 0.5 at r ~ a 
few r_2; then (3 decreases outwards, though with wide 
oscillations (e.g., Navarro et al. 2010; Ludlow et al. 
2010). Remarkably, within r_2 the anisotropy param- 
eter and the slope of the density profile appear to be 
correlated (sec Huss, Jain & Steinmetz 1999; Hansen & 
Moore 2006) through a /3 — 7 relation of approximate 
form /3(r) « -0.15 -t- 0.27(r). 

Finally, DM halos show a powerlaw spherically- 
averaged profile of the quantity Q{r) = p/cr^. This has 
the dimensions of a phase-space density, but it is not a 
true measure neither of it nor of its coarse-grained ver- 
sion (see discussion by Ascasibar & Binney 2005; Sharma 
& Steinmetz 2006); thus it is often referred to as a 
pseudo (or a proxy of the) phase-space density, although 
it is still debated whether the radial or the total veloc- 
ity dispersion should enter its definition (see discussion 
by Schmidt, Hansen & Maccio 2008). Remarkably, al- 
though the density p{r) and velocity dispersion iJ^(r) 
have articulated runs, Q(r) follows a simple powerlaw 
Q{r) oc over three order of magnitude in radius, 
with the same exponent x « 1.9 applying for all halos 
(Taylor & Navarro 2001; Hoffman ct al. 2007; Ascasi- 
bar & Gottlober 2008; Vass ct al. 2009; Navarro ct al. 
2010; Ludlow et al. 2010). On the other hand, the re- 
cent simulations highlight that such a powerlaw behavior 
holds within r_2, but that Q{r) possibly steepens in the 
inner regions, while it flattens appreciably in the out- 
skirts. The origin of the powerlaw behavior is presently 
unknown, but numerical experiments with different per- 
turbation spectra indicate that it is not related to initial 
conditions or hierarchical merging (see Wang & White 
2009). 

Unfortunately, these remarkable flndings still lack a 



firm theoretical background, while their reliability is re- 
stricted in the inner and outer regions of the halos by lim- 
ited resolution and small particle statistics, respectively. 
These are good reasons to complement the numerical ap- 
proach with analytic models. 

In the analytic vein, a simple approach is provided 
by the self-gravitating, static equilibria of DM based on 
the Jeans equation (Taylor & Navarro 2001; Austin et 
al. 2005; Dehnen & McLaughlin 2005; Lapi & CavaUere 
2009a, 2009b, 2011) 

Id. ^^(Jr GM(< r) „ 

p dr ^ ^2 

Two ingredients are needed to solve it for p{r): the profile 
of the anisotropy parameter /3(r), and that of the phase- 
space density Q{r) linking p(r) and cr^(r) in the spirit 
of an 'equation of state'. Both runs can be predicted 
from simple scaling arguments, and refined by compari- 
son with numerical simulations (sec above) , to the effect 
that l3(r) takes the form of the linear (3 — ^ relation, 
while Q{r) oc r~^-^ applies. Interestingly, the resulting 
density profiles are well described by the Sersic-Einasto 
formula (Lapi & Cavaliere 2011), and turn out to be very 
close to the simulation results in the halo middle where 
the Jeans solutions are reliable; conversely, Ludlow et 
al. (2011) take up from simulations an Einasto shape for 
p{r) and solve Jeans to closely recover the powerlaw run 
of Q{r). However, such analytic approaches are limited 
being based on a static equation; they provide final equi- 
librium pictures of halos, but little can tell on how DM 
particles progressively collapse or accrete, and then re- 
lax. As such, they sidestep origins and building up of the 
profiles Q{r) and /3(r), and arc expected to fail close to 
the halo outskirts, in the region exposed to infall where 
equilibrium is not yet achieved. 

To go beyond these limitations requires formulating 
physical models of infall and accretion in an expand- 
ing universe; to that purpose we resort to an advanced 
fluid-like description of the DM dynamics including ran- 
dom non-radial motions and relaxation effects, but still 
amenable to a self-similar treatment. The latter pro- 
vides a tractable, analytic way of investigating time- 
dependent problems in complex physical systems; it is 
possible whenever the system dynamics can be charac- 
terized, besides spacetime variables, by a few parameters 
with independent dimensions. Although providing only 
a particular analytic solution to the physical problem, it 
often accurately yields the long-time behaviors and of- 
fers a useful guide for understanding the generic features 
of the system (see discussions by Sedov 1959; Zeldovich 
& Raizer 1967). In the present context, a self-similar 
description for the collapse and condensation of a DM 
halo is allowed by choosing a scale-invariant shape of the 
initial mass perturbation, and adopting an Einstein-de- 
Sitter cosmological framework. 

Self-similar solutions for the purely radial infall of col- 
lisionless matter have been pioneered by Gunn & Gott 
(1972), then explicitly derived by Fillmore & Goldreich 
(1984) and Bertschinger (1985) on adopting a Lagrangian 
description for the orbits of individual particles, and by 
Teyssier et al. (1997) and Subramanian (2000) on adopt- 
ing an equivalent fluid approximation for the ensemble of 
particles. Such self-similar treatments with purely radial 
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infall concur in yielding very steep inner density profiles 
with 7 '-^ 2, at variance with the simulation outcomes. 

White & Zaritsky (1991), Ryden (1993), Sikivie et al. 
(1997), Subramanian (2000), Nusser (2001), Zukin & 
Bertschinger (2010a, 2010b), Vogelsberger et al. (2011), 
and Lithwick & Dalai (2011) took steps toward fixing the 
problem by inclusion of non-radial motions originated 
by tidal torques, either in spherical or triaxial collapses; 
when angular momentum is present, the system develops 
a tangential velocity dispersion that causes a flattening 
of the inner density profile (see also Eq. 1). However, 
the quantitative effect is sensitive to the specific amount 
of angular momentum endowed at, or acquired during 
the infall, so the solutions still fail to explain the ap- 
proximately universal shapes found in simulations (see 
discussion by Lu et al. 2006). 

The view we submit here envisages that all such short- 
comings go back to a key ingredient missing (or at least, 
inadequately treated) in the classic self-similar mod- 
els, i.e., dynamical relaxation. This involves two differ- 
ent mechanisms leading to broadly similar macroscopic 
outcomes: phase-mixing, a localized process related to 
spreading of neighboring particle orbits over phase-space; 
and violent relaxation, a volume process driven by ir- 
regular fluctuations of the gravitational potential. Both 
these effects have long been recognized to play a cru- 
cial role in the approach of a gravitational system to its 
equilibrium configuration, with violent relaxation espe- 
cially relevant when the initial condition of the collapse is 
clumpy (e.g., van Albada 1982; Ma & Bertschinger 2004; 
Lu et al. 2006). This is certainly the case in the stan- 
dard cosmological framework, where DM objects form 
from accretion punctuated by mergers. 

For example, Lynden-Bell (1967) and Nakamura 
(2000) based on a statistical mechanical approach to 
show that an equilibrium state may be achieved when 
violent, full relaxation erases the memory of the initial 
conditions. Besides delicate issues concerning the math- 
ematical consistency of these theories (see discussion by 
Arad & Lynden-Bell 2005), a number of numerical exper- 
iments have demonstrated that the relaxation is always 
incomplete (see van Albada 1982; Henriksen & Widrow 
1999; Trenti, Bertin & van Albada 2005, and references 
therein), with a significant correlation retained between 
the final and initial state of a particle. In fact, this may 
explain why the detailed density profiles of virialized DM 
halos ultimately do depend on the formation history and 
on cosmogonical conditions. 

In this paper we propose an advanced treatment of 
both dynamical relaxation effects within self-similar solu- 
tions for the spherical collapse of DM halos in a fluid-like 
approach. We will show that our description resolves the 
discrepancy between the outcomes of self-similar models 
and numerical simulations, in that the solutions of the 
former feature profiles of mass density, pseudo phase- 
space density and anisotropy consistent with the latter. 

The plan of the paper is straightforward. In § 2 we 
derive the evolution equations for the moments of the 
DM distribution function in the fluid approximation, in- 
cluding the terms responsible for dynamical relaxation 
(see Appendix). In § 3 we reduce these equations to a 
self-similar form. In § 4 we derive conservation laws, 
investigate the inner asymptotic behaviors of the solu- 
tions, present the full solutions, and compare them with 



the simulation outcomes. Finally, in § 5 we discuss our 
approach and summarize our main conclusions. 

2. DM DYNAMICS IN THE FLUID APPROXIMATION 

We consider a spherically symmetric, nonrotating DM 
halo, and describe particle positions and velocities in 
terms of the standard polar coordinates r, 9, (j), and 
of the corresponding velocity components Vr, ve, v^. 
In the phase-space the particle distribution function 
f{r,Vr,V0,v^) evolves as dictated by the classic Boltz- 
mann equatiorQ 

dt f+Vr dr f+Vr dy^ f+ve dyg f+v^ d^^ f = {dt /)* . (2) 

Here {dt /)* is a term that describes collective 'collisions' 
contributing to dynamical relaxation, and will be spec- 
ified in § 2.1 and in the Appendix. The single-particle 
velocity components change according to the equations 
of motion 

GM(<r) , v^g+vl 



-- 



Vr V0 



Vr Vcjy 



r tan( 



Ve Vcf, 
r tanfl 



(3) 



where M{< r) is the DM mass within a radius r. 

Next we derive the evolution equations for moments of 
the velocity distribution, up to the second order inclusive. 
In the following the 'mean' value < X > of a generic 
quantity X is defined by 



<X>- 



1 



<rvfX , 



(4) 



where the spatial density is p = Jd'^vf. In the tan- 
gential directions, the mean components of the velocity 

< Ve ^ij, >= must vanish (as all other moments of odd 
order) while the corresponding dispersions = cr^ must 
be equal (in fact, the value of any moment is unaffected 
by interchanging ve and w^). In the radial direction, in- 
stead, one has to consider the mean velocity u =< Vr > 
and the dispersion cr^ =< {vr — u)"^ >=< v'^ > —v}. 

To obtain a closed set of moment equations we need 
to relate the third moment with the lowest-order ones. 
The standard and simplest assumption is to require 
zero radial skewness < (vr — uf' >— 0, which yields 

< vf >= u{v? -I- 3tT^); in the literature this is often 
referred to as ^ fluid approximation' (see Teyssier et al. 
1997; Chuzhoy & Nusser 2000; Subramanian 2000; also 
Dehnen & Read 2011 and references therein). Radial 
skewness, on the other hand, would imply an outward 
flow of energy through the system, or in other words a 
tendency for the most energetic particles to move prefer- 
entially outwards. Skewness is often taken into account 

^ For the sake of simplicity, in spherical symmetry and in the 
absence of net macroscopic rotation of the halo, one can neglect 
any explicit dependence of / on the angular variables d and </> (see 
Larson 1969, 1970; Subramanian 2000). In fact, the related terms 
{dgf)/i", {d(f,f)/r sin6 to be added on the l.h.s. of Eq. (2) would 
introduce in the moment equations (Eqs. 5 to be derived next) 
averages like < vg ^ > or dispersions like a^g that must 

vanish anyway in a macroscopically spherical, nonrotating halo. 
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in studies of star clusters' outskirts to evaluate the sec- 
ular evolution of the system due to the escape of ener- 
getic stars; on the other hand, two-body encounters tend 
to erase skewness at the center over a timescale slightly 
longer than the local relaxation time, see Lynden-Bell & 
Wood (1968) and Larson (1970). These analyses suggest 
skewness to be smaller than second-order moments by 
one order of magnitude or more throughout the system. 
In the context of DM halos accreting from a cosmologi- 
cal environment, it has been shown (see references at the 
beginning of the paragraph) that skewness is negligible 
in the inner region where large amount of phase-mixing 
has occurred, and it may become relevant only in a nar- 
row region aroimd the virial boundary (or better, the 
outermost caustics) , which in the fluid approximation is 
considered as a sharp discontinuity and dealt with as dis- 
cussed in § 3.1. 

Multiplying Eq. (2) by 1, Vr, {vr ~ u)^, Vg and inte- 
grating over the velocity distribution yields the system 
of coupled partial differential equations 

{dt+udr)p+-^dr{r^u)=0 , 

{dt+udr)cjl+2aldrU={dt(J% , (5) 

{dt +udr)a0 + '^agu = {dt ci^)* , 
drM = 47rr^p . 

These are commonly referred to as the continuity, mo- 
mentum, energy, angular momentum, and mass equa- 
tions; note that the second constitutes the extension of 
Eq. (1) to time-dependent conditions. 

Remarkably, the first and the third Eqs. (5) can be 
combined into the form 

id, + udr)l0g^=^^^. (6) 

This describes an entropy flow for a one-dimensional fluid 
with effective density p = pr^, effective pressure p = 
pcF^i and effective specific entropy s = \ogp/p^, in terms 
of a microscopic adiabatic index equal to 3. 

2.1. The 'collision' terms 

We now turn to discuss the terms (9t crj: g)j, on the 
r.h.s. of the third and fourth Eqs. (5), that describe col- 
lective 'collisions' in the DM fluid; hereafter these will 

be simply referred to as collision terms. Physically, they 
arise because dumpiness in the infalling matter induces 
fluctuations of the gravitational potential that collec- 
tively affect the dynamics of DM particles. In fact, a 
particle effectively sees a stochastically fluctuating poten- 
tial, and (consistent with the fluctuation-dissipation the- 
orem) suffers of a dissipative drag that produces an effec- 
tive 'cooling' (see Chandrasekhar 1943; Kandrup 1980; 
Antonuccio-Dclogu & Atrio-Barandela 1992; Del Popolo 
& Gambera 1997; Ma & Bertschinger 2004; for a re- 
view see Binney & Tremaine 2008). On the other hand, 
as pointed out by Bekenstein & Maoz (1992) and Ma & 



Boylan-Kolchin (2004), effective 'heating' of a particle by 
the fluctuations can also constitute a relevant process. 

In the present context, we consider a halo of mass M 
including clumps of typical mass m. In Appendix A we 
work out the form of the collision terms in the Fokker- 
Planck approximation, under the assumption of a quasi- 
isotropic distribution function; the results are as follows 

{dta^^)^ = f— G^mplogA " 3 ^ , 

(7) 

{dt al). = +^a^ G'mp log A , 

with log A = logAf/m the standard 'Coulomb log- 
arithm' (see Appendix A. 2 for details). Note that 
the ratio of the average relaxation timescale ~ 
10"^ cr^/G^ mp log(M/m) from Eqs. (7) to the Hubble 
time (3/87rGp)^/^ reads ~ a few lO'^A/"/ \ogM in terms 
of the effective clump number J\f = M/m; this implies 
that relaxation can be efficient for J\f < 10, i.e., when 
a limited number of clumps is present (a more detailed 
discussion of the relaxation strength is provided in Ap- 
pendix A. 2 and recalled in § 3). 

Note that in Eqs. (7) the (negative) cooling term domi- 
nates in the radial direction, while the (positive) heating 
term dominates in the tangential direction. The colli- 
sions induced by clumps sharing the accretion inflow do 
not change the total velocity dispersion a^^^ = cr^ -|- 2 fjg, 
since (di (Jj^ot)* = holds; this means that the overall ran- 
dom component of the energy is conserved, being just re- 
distributed between the radial and the tangential degrees 
of freedom. On the other hand, collisions tend to erase 
any velocity anisotropy since [dt {af — cr^ )]* oc — {af — aj) 
holds. 

Such collective collisions provide one mechanism con- 
tributing to dynamical relaxation, close to the classic but 
still incomplete notion of violent relaxation (see § 1 and 
references therein). The other mechanism is related to 
phase-mixing, and its implications in the present context 
will be discussed in § 3.1. 

3. SELF-SIMILAR DESCRIPTION 

Self-similar solutions of the above Eqs. (5) obtain un- 
der two assumptions: (i) an Einstein-de-Sitter cosmo- 
logical framework, which still provides an increasingly 
good approximation to the concordance cosmology for 
z > 0.5 where most galactic halos form and evolve; (ii) a 
spherically-symmetric, scale-free shape of the initial DM 
mass perturbation of the form 

(8) 

M 

which may be considered as a piecewise approximation 
to a realistically bell-shaped cold DM perturbation. Here 

SM rciprcscnts the mass excess in a shell of initial co- 
moving radius ri oc M^/^ enclosing a mass M of mat- 
ter at background density. Such a shell will progres- 
sively detach from the Hubble flow, reach a maximum 
'turnaroimd' radius i?ta oc ri/{5M/M) oc iW^+^Z^, and 
collapse back to the standard virial radius i?200 ^ -Rta/2. 

The virialization occurs when 5M/M attains the crit- 
ical threshold 1.686 Z)~^(t) in terms of the linear growth 
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factor D{t) oc i^/^ along the cosmic time t in the 
Einstein-de-Sitter cosmology. So the shape parameter e 
also governs the mass buildup after M{t) oc D^^'^{t) oc 
t^/^^, or equivalently sets the collapse time tcoii = 
M/M = 3et/2 for the shell surrounding the mass M. 
Values e < 2/3 correspond to a stage of fast collapse 
relative to the Hubble expansion, while values e > 2/3 
correspond to a slow mass accretion. 

The scaling M cx t^^^'^ oc (1 + z)~^^'^ may be conve- 
niently compared with that of the characteristic cluster- 
ing mass Mc oc (1 + 2)-6/(n+3) scale-free hierarchi- 
cal cosmology with a power spectrum P{k) oc A;" (e.g., 
Peebles 1980). Perturbations characterized by a specific 
value of e therefore accrete mass at the same rate as 
the 'typical' mass with e = (n + 3)/6. For example, 
the fast collapse of the perturbation bulk for a galactic 
halo is described bye« 1/12 -^1/6 corresponding to 
nRi —2.5-=-— 2, while for the halo of a galaxy cluster val- 
ues e « 1/6-^1/3 apply corresponding to n w — 2 -= — 1. 
On the other hand, during the development of the out- 
skirts, the accretion rate slows down and e can grow con- 
siderably larger (see discussion in § 5). 

Once the shape of the initial perturbation has been 
specified, there are no additional physical scales in the 
problem, and the timc^ evolution of the system must ap- 
proach self-siniilarity after a short initial transient. This 
implies that a single solution describes the structure and 
time behavior of the system, when expressed in properly 
scaled variables. Since the current turnaround radius is 
easily seen to grow as 



Rta.it) oc , with 



C-^(l + 3e) 



(9) 



it is convenient to introduce the self-similar variable 
A = r/Rtait), and define the adimensional radial veloc- 
ity, density, mass, and velocity dispersions through the 
relations 



u(r, t) 



Rk 



U{X) 



p{r,t) = 



V{X) 



(10) 



Thus we can transform Eqs. (5) into the system of 
ordinary differential equations 

1 



(u-xou' + {^-i)u + ^{vi:iy + 

+ ^S2-E2)+^^-0 



(11) 



{u-xoi^iy+2 (?-i + 

M' = 3X'^V , 



U 



where prime denotes differentiation with respect to A. 

Moreover, with a constant value for the effective num- 
ber Af = M/m of clumps in the infalling matter, the 
collision terms are seen to scale as 1/t, and hence they 
do not break the self-similarity. Using the variables de- 
fined in Eqs. (10) we can write them in the form 



-'tot 



(12) 



ready to be inserted in Eqs. (11). In the above expression 



(13) 



QAVStt log J\f 



405 



is the strength parameter of the collective collisions, 
which depends mainly on the 'dumpiness' 1/J\f of the 
infalling matter, see Appendix A and Fig. Al for details. 

In solving numerically the fully self-similar Eqs. (11), 
we will consider different values of k^, guided by the fol- 
lowing physical considerations. During the fast accre- 
tion (e < 2/3) we expect that a limited number of major 
clumps Af <10 rapidly merge to build up the halo bulk; 
after Eq. (13) this implies efficient dynamical relaxation, 
with strength parameter taking on values k^, w 0.1. On 
the other hand, during the late development of the out- 
skirts we expect slow accretion (e > 2/3) of many small 
clumps with Af > 30; these conditions imply inefficient 
dynamical relaxation with < 0.01. 

We stress that such values for the effective number of 
clumps in the different accretion regimes are consistent 
with the findings of numerical simulations (see Wang et 
al. 2011, their Fig. 7). The latter show that most of the 
mass in the inner region (which in turn is about 20% of 
the total) is contributed by a number of 5 — 10 major 
mergers with mass ratios 1 : 2 — 1 : 3; on the other hand, 
most of the mass in the outskirts is contributed by a 
number > 30 of minor mergers with mass ratio < 1 : 10, 
and the rest by smooth and diffuse accretion. 

3.1. Boundary conditions 

Since the original Eqs. (5) describing the DM flow are 
hyperbolic, discontinuities are expected to develop. How- 
ever, given the collisionless nature of the DM particles (at 
the binary level), the discontinuities are constituted not 
by shocks, but rather by caustics where the bulk infall en- 
ergy is partly converted into random motions. The caus- 
tics occur very close to the turning points where the ra- 
dial velocity of a shell vanishes (e.g., Bertschinger 1985); 
there adjacent mass shells catch up with each other, and 
the intervening matter is compressed to a divergent den- 
sity. Inward of a caustic multiple 'shell crossings' occur, 
the effective gravitational force experienced by a particle 
is (albeit slowly) time dependent, and neighboring parti- 
cles' orbits go out of phase, causing the so called 'phase- 
mixing'. Deeper in radius and/or later in time, more 
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shell crossings have occurred, to originate a smoother 
coarse-grained particle distribution in phase-space. 

In the treatments based on the Lagrangian viewpoint 
where the orbits of single DM particles are followed, 
phase-mixing primarily takes place in an outer layer in- 
cluding a few caustics (see Fillmore & Goldreich 1984, 
Bertschingcr 1985). In the fluid approximation one ren- 
ders the layer of such caustics with a single "discontinuity^ 
across which the relevant Rankinc-Hugoniot— type jump 
conditions apply. An extensive literature has shown the 
effectiveness of such an approach in closely matching the 
single-particle results (see Teyssier et al. 1997; Chuzhoy 
& Nusser 2000; Subramanian 2000). On the other hand, 
in the volume inward of the discontinuity the other re- 
laxation mechanism related to collective collisions can 
operate (see § 2.1). 

The radius where the caustic discontinuity occurs in 
the self-similar description must be a constant fraction 
Ac of the turnaround radius, and will be computed below 
as an eigenvalue (see Bcrtschinger 1985). In detail, we 
proceed as follows. First of all, outside Ac the evolution is 
identical to the turnaround and collapse of a pressureless 
mass shell; a parametric form of the upstream accretion 
flow is given by (e.g., Peebles 1980) 



A = sin^ 



e /e-sinl 



U{\) 



A sin 6* (0- sin 6*) 
(1- cos 61)2 ' 

(61 -sin 61)2 



(14) 



A^(A) 



2 (l-cos6»)3 (l-h3ex) 

9 3 (61 -sin 0)2 
2 (l-cos6')3 ■ 



3 U(X) 

with x=l---^ 



Note that at the turnaround position corresponding to 
A = 1 the above Eqs. (14) give U{1) = S^_g(l) = 0, 

M{1) = (37r/4)^ and V{1) = Al(l)/(1 + 3e). ' 

Then, by integrating the self-similar Eqs. (11) in a 
region across the caustic discontinuity, we obtain the 
jumps 



1 



{U- - CAe) 



also generated by a nonsphcrical collapse, but a simple 
treatment like the present one cannot tell to what degree 
this occurs; in fact, no source term for the tangential 
dispersion is present in the corresponding moment equa- 
tion, meaning that it must be assigned as a boundary 
condition. For example, Zukin & Bertschingcr (2010a) 
took steps toward estimating the tangential dispersion 
at the turnaround basing on tidal torque theory; their 
result overestimates somewhat (by a factor about 2) the 
outcomes of iV— body simulations, and ultimately must 
be tuned in terms of the latter. This issue is beyond 
the scope of the present paper, and we just parameterize 
the tangential dispersion at the caustic discontinuity in 
terms of the boundary value /3+ = 1 — (Eg/E2)^ for the 
anisotropy parameter. We shall find that such a value 
does not affect the inner shapes of the mass and pseudo 
phase-space density profiles when efficient dynamical re- 
laxation is at work. 

Technically, the location of the caustic discontinuity Ac 
is determined as an eigenvalue (e.g., Bertschingcr 1985), 
by imposing the following inner physical constraints 



A^(0) =ZY(0) =0 , 



(16) 



i.e., the mass and velocity at the center must vanish. 
Note that the caustic constitutes an effective boundary 
for the halo, which is close if exterior to the virial radius 

-^200- 

4. SELF-SIMILAR SOLUTIONS 

We now turn to solving the self-similar Eqs. (11). Be- 
fore handling the problem numerically, useful insights 
into the behavior of the solutions is found from the ana- 
lytic work that follows. 

4.1. Conservation laws 

First of all, we find that integrals of motion are asso- 
ciated to the self-similar Eqs. (11) when collisions terms 
are neglected. It is convenient to introduce the auxiliary 
function (see Chuzhoy & Nusser 2000) 



= exp 



dC' 



(17) 



V, 



2V. 



(15) 



so that U = + TIT' and ZY' = C + 1 - TT"/{T'f 
obtain. Then from the continuity, mass, energy, and an- 
gular momentum equations we find 

Voi.X-'^T^-^^T' , 



16 



here the — , + signs refer to upstream and downstream 
values, respectively. Outward of the discontinuity, we 
have taken the radial and tangential velocity dispersion 
to be null, since only bulk radial motions are associated 
to the inflow of a shell until it collapses back to Ac and 
undergoes shell crossing. At the discontinuity, the flow 
is compressed and the bulk radial velocity is converted 
into radial dispersion, similarly to what occurs with the 
usual jump conditions at a shock but for a fluid with one 
degree of freedom. 
On the other hand, tangential velocity dispersion is 



M oc 



2-3^ 



T- 



•2-3 C 



oc T-^^ {T'f 

Hlo^X-^T^^^-^^^ . 
Rearranging these relations yields 



(18) 



M oc 



2-3^ 



X^V{U- AO , 



e: 



A4 P2 



^ OC , 



(19) 
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which constitute the 'conservation' laws of mass, effective 
entropy, and angular momentum. In fact, like in Eqs. (6), 
the quantity logp^/A* I?2] ^ \og[(V T.^) / (V X^)^] 
may be interpreted as an effective entropy for a one- 
dimensional fluid with effective density VX^, effective 
pressure (2?A^)E^, and microscopic adiabatic index 
equal to 3. Since the mass is a monotonically increasing 
function of the radius, it follows that when e < 2/3 ap- 
plies such an effective entropy increases toward smaller r; 
physically, this is interpreted as an efficient relaxation of 
the particles during the fast collapse of the halo inner re- 
gions. On the other hand, when e > 2/3 applies the effec- 
tive entropy grows with r; physically, this is interpreted 
as a progressive stratification of the particles' orbits (or 
better, of the orbit apocenters) during the slow accretion 
that builds up the halo outskirts (see Bertschinger 1985; 
also Taylor & Navarro 2001 for a similar discussion in 
terms of the pseudo phase-space density). 

When the collision terms are efficient, only the mass 
conservation still applies throughout the halo. 

4.2. Asymptotic behaviors 

We now derive analytically the asymptotic behavior 
of the solutions near A ~ 0. For the sake of simplicity 
and with no loss of generality, we assume the following 
powerlaw forms of the mean radial velocity, density, and 
velocity dispersion^ 



U ^UqX 



Vr^VnX- 



(20) 



in terms of two exponents 7 and w; correspondingly, the 
mass behaves as ~ MoX^~^ withA^o = 3I?o/(3-7), 
and the pseudo phase-space density Q = T) jY? follows 
Q ~ Qo with So = 2?o/S3 and X = 7 + 3a;/2. 

For such exponents, the continuity and energy equa- 
tions yield the relations 



7 



-2-1-3^0 



(21) 



w = -2 



In addition, when collision terms are neglected, the mo- 
mentum equation writes 



AWo (Z^o-l) + (S?)o A"-i (-7+c.+2/3o) + ^ ^ 



7 



A^-'' -0: 
(22) 



^ Given the constraint W(0) = from Eq. (16), in principle one 
should write W ~ A" with > 0. However, it can be shown that 
Eqs. (11) do not admit asymptotic solutions for v < \\ hence for 
u > 1 one can write W ~ a-t the first order. Moreover, for u > 1 
the energy and angular momentum equations imply the radial and 
tangential dispersions to scale asymptotically in the same manner, 
so that one can write „ ~ A". 



hereafter /3o stands for the central value of the anisotropy 
parameter. 

There are two possibilities for satisfying Eqs. (21) and 
(22). 

• First, in Eq. (22) the exponents of the middle and last 
terms are negative and equal, while the coefficient of the 
middle term must be negative; then uj = —7 + 2 applies, 
and 7 > 1 + (3q must hold. Together with Eqs. (21), these 
yield VIq ~ and 



7 



9e 



l-|-3e 



2 _ 2 - 3e 
1^1+3^ ' 



(23) 



X 



3 2 + 3e 
2 1 +3e 



for e > 



1 l + /9o 
3 2-/3o 



This case corresponds to the late regime of slow accretion 
of DM particles onto a preformed halo bulk; in fact, the 



inner density in physical units p(t) oc X t 



i-(2-7«) 



behaves as pit) cx i° so that it is independent of the time 
t, i.e., is unaffected by the outskirts growth, consistent 
with the two-stage formation picture. 
• Second, in Eq. (22) the coefficient of the middle term is 
zero so that — 7-|-a;-|-2/3o = applies, while the difference 
2— 7— w > between the exponents of the last and middle 
term must be positive to imply 7 < l-|-/3o- Together with 
Eqs. (21), these relations yield Z^o - 4 [1 - 6e + /3o (1 + 
3e)]/[9e(2/3o-5)] and 



7 



X 



3e + 2 + 2^ 
^ 3e + 7 ' 

3(3e + 2)-2/3o (Se -I- 4) 
(3e + 7) 

3 (3e + 2) (5-2;9o) 
2 3e + 7 



(24) 



for e < 



1 l + /3o 
3 2-/3o 



This case corresponds to the early fast collapse of the 
halo bulk; in fact, the inner density in physical units 
behaves as p{t) oc X'^ t'^ oc t4 [i-6.+/3o (3e+i)]/3e (3.+7) 
so that it grows with the time t. We remark that 
the asymptotic relation 7 = w + 2 /?o holds like in the 
solutions of the anisotropic Jeans equation found by 
Dehnen & McLaughlin (2005) and Lapi & Cavaliere 
(2009b). We also stress that the expressions of the 
asymptotic slopes in Eqs. (23) and (24) are continuous 
at e = (l+/3o)/3(2- /3o). 

Note that the classic results based on assuming adi- 
abatic invariance of the radial action may be recovered 
in the fluid approximation by imposing the solution to 
feature for A ~ the additional regularity condition 
tdrU — U'{X) ~ 0, or equivalently ~ Q (see Teyssier 
et al. 1997, Subramanian 2000). While the slope 7 given 
by Eq. (23) for e>{l+ ,9o)/3 (2 - /3o) is consistent with 
~ 0, that given by Eq. (24) for e<{l + /3o)/3 (2 - /3o) 
is not as Z^o < holds; thus for any e in the latter range, 
the maximal slope 7 = 1 -I- /3o consistent with the condi- 
tion Z//o — is to apply. We give three relevant examples: 
in a purely radial collapse with /3(r) — j3a = 1 the condi- 
tion Wo — would imply 7 = 9 e/(l -I- 3e) for e > 2/3 and 
7 = 2 for any e < 2/3, which is the Fillmore & Goldreich 
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(1984) result; in an isotropic core with /3o = it would 
imply 7 = 9e/(l + 3e) for e > 1/6 and 7 = 1 for any 
e < 1/6; for a slope 7 < 1 to hold at the center it would 
imply /3o < 0, i.e., prevailing tangential motions (Subra- 
manian 2000). We remark that although not satisfying 
the condition Z//o — 0, our asymptotic solution given by 
Eqs. (24) is physical and suited to describe the fast col- 
lapse of the halo bulk; to the best of our knowledge, this 
was not previously known. 

In Fig. 1 we plot against e the inner asymptotic be- 
haviors given by Eqs. (23) and (24) for the slopes of the 
density 7 and pseudo phase-space density x, for different 
values of /3o. We remark that the full solution attains 
its asymptotic, shape quite slowly, as a result of a log- 
arithmic convergence; e.g., the solution with e = 1/6, 
that features a central density slope —1, has still a slope 
around —1.3 at r w 0.1 r_2- 

When collision terms are efficient, Eqs. (23) and (24) 
are still valid, but the asymptotic behavior 



e ~ exp 



■IAq (5]tot)Q 



A' 



3-27-3w/2 



47-t- 3w - 6 ^ 

(25) 

applies, and enforces /3o = 1 — (S^)o/ (J^Do — to hold at 
the center. Then from Eqs. (23) and (24) the value e k, 
1 /() is seen to sciparatc; the violent collapse of the inner 
region from the calm, inside-out growth of the outskirts. 
We stress that dynamical relaxation in the inner region 
is mainly provided by collective collisions, whilst in the 
outskirts it is related to phase-mixing (see also § 5). 

Note that cfhcicnt dynamical relaxation during the 
fast collapse stage enforces a vanishing central anisotropy 
^0 = 0, while making the inner mass and pseudo phase- 
space density profiles only weakly dependent on both the 
perturbation shape parameter e < 1/6, and the outer 
anisotropy parameter In this sense, the inner halo 
structure turns out to be approximately universal. 

4.3. Numerical solutions 

We solve numerically the system of ordinary differ- 
ential Eqs. (11) over the spatial range 10"'* < A < 1 
with an Adams-Bashford-Moulton method of variable or- 
der, adaptive stepsize, and error control; the location of 
the caustic discontinuity is determined with a standard 
shoot-and-match technique, by requiring the solution to 
satisfy the inner constraints Eqs. (16), while the jump 
conditions Eqs. (15) are applied across the discontinuity. 

As a preliminary check, in Fig. 2 we present the solu- 
tion for e = 1, /3+ = 1, and no collision terms (k* = 0) 
corresponding to pure radial infall onto a point-mass per- 
turbation; this is the classic case solved by Fillmore & 
Goldreich (1984) and Bertschinger (1985) basing on self- 
similar particle trajectories, and cquivalently by Tcyssier 
et al. (1997) and Subramanian (2000) in the fluid ap- 
proximation. We recover in detail the solutions of the 
latter authors. 

Then we focus on the vahic c = 1/6 that corresponds 
to a spectral index n = —2 typical of a galactic halo. 
Figs. 3-4-5 refer to /3+ = 1, /3+ = 0.5 and /3+ = 0.25, re- 
spectively, still in the absence of collisions (k* = 0); these 
illustrate how the inner density profile is fiattencd rela- 
tive to the purely radial case when non-radial motions 
are included. However, as already stressed by Subrama- 
nian (2000), Nusser (2001), and Zukin & Bertschinger 



(2010b) such a flattening depends on the amount of an- 
gular momentum assigned at the caustic discontinuity 
and on how mass shells are torqued after turnaround, 
so that producing the approximately universal shape of 
the inner density profiles found in simulations requires 
fine tuning of a sort in the initial conditions and/or in 
the torque mechanism. Moreover, note that in the ab- 
sence of collision terms, the anisotropy profile /3(r) rises 
inward from the boundary value a behavior at strong 
variance with what is seen in numerical simulations. 

In Figs. 6-7-8 we retain the values e = 1/6 and j3+ = 
0.25, but include efficient collision terms; from Fig. 6 to 
7 to 8 the strength parameter k^, of collisions is increased 
from 0.01 to 0.05 to 0.1. The effect of coUisions is twofold: 
the inner slope of the density profile is now flattened to 
values 7 < 1, while the anisotropy parameter /3(r) is 
lowered inward to a vanishing value j8o ~ 0. As collisions 
become more and more efficient, these effects occur on 
wider and wider scales. The same qualitative bcihavior 
takes place for other values of e < 1/6, as illustrated 
in Fig. 9 for the specific case e = 1/8, and still with 
/3+ = 0.25 and = 0.1. 

In Fig. 10 we plot the position of the caustic disconti- 
nuity Ac as a function of e, for two different values of 
with and without efficient collision terms. At fixed /3+ 
and it is seen that Ac increases with e; this is because 
as e grows and the infall rate slows down, the lower in- 
fall stress allows the caustic discontinuity to be located 
farther out. At fixed e and k*, a lower value of corre- 
sponding to a higher angular momentum, yields a larger 
value of Ac; this is because for particles with higher an- 
gular momentum it is harder to penetrate deep into the 
halo. Finally, at fixed e and /3+, a higher k^, yields a 
smaller Ac since collisions isotropize the velocity disper- 
sions, and so the particles lose part of their initial angular 
momentum and can penetrate deeper into the halo. 

We stress that all the above self-similar solutions fea- 
ture a wide radial range from the center to about a few 
r_2 where the bulk velocity u is approximately null. 
This implies that in the second of Eqs. (5) the term 
{dt -h u dr) u ~ closely vanishes, and the equation itself 
reduces to a Jeans form like Eq. (1). As such it describes 
a nearly static equilibrium, endowed with runs of Q{r) 
and /3(r) as provided by the full system of equations. 

4.4. Comparison with N— body simulations 

In Figs. 11-14 we compare our self-similar solutions to 

the outcomes of state-of-the-art numerical simulations. 
In all these plots, for the solution inward of r_2 our fidu- 
cial values arc e = 1/8 (as representative for the range 
e < 1/6) and = 0.1; in fact, these apply to the fast in- 
ner collapse of galactic halos with effective spectral index 
n = 6e — 3 « —2.25, when rapid merging of an effective 
number of clumps A/" ^ 10 implies efficient dynamical re- 
laxation (see § 3 and Fig. Al). For the solution outward 
of r_2, our fiducial values are e = 1/2 (as representative 
for the range e > 1/6) and = 0.01. This is because we 
expect the outskirts growth to be dominated by smooth 
accretion from the tapering perturbation wings with ef- 
fective spectral index n = 6e — 3 « 0; now the accre- 
tion involves many small clumps with A/" > 30, implying 
dynamical relaxation to become inefficient. We stress 
that similar behaviors obtain in the two radial ranges for 
reasonable variations of e and around our reference 
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values. The overall picture to be compared with real or 
simulated data may be obtained on matching these two 
solutions around r_2. 

As to the other parameter this is less amenable 
to physical pinpointing. We expect /?+ < 1, i.e., a de- 
viation from purely-radial collapse, due to asymmetries 
both in velocity and configuration space. On the other 
hand, from the experimentations reported in § 4.3 we 
know that in the presence of k-^ > 0.01, specific bound- 
ary values of /?+ do not materially affect the solution, 
including the run of /3{r) vanishing toward the center. 
Finally, note that around the outer caustic discontinu- 
ity /3 = 1 — Eg /E^ is ill defined anyway, since both the 
radial and the tangential dispersions are small upstream 
(see also end of this section). Within the above con- 
straints, we adopt the value /3+ = 0.25 both for the inner 
and outer solutions, because around r_2 this yields sim- 
ilar values of (3, and provides a close match of the two 
(already comparable) density slopes. 

In Fig. 11 we compare the inner self-similar profiles of 
density (top panels) and density slope (bottom panels) 
with standard fitting formulae to the outcomes of nu- 
merical simulations. Specifically, we illustrate the NFW 
profile (dot-dashed line) and the Einasto profile with 
shape parameter rj = 0.17 (dashed line); in the inner 
radial range r > 0.04 r_2 that is currently accessible to 
numerical simulations, the former constitutes a popular 
description of the virtual data, while the latter has been 
recently found to constitute a better functional represen- 
tation (see § 1). It is seen that the inner self-similar so- 
lution follows closely the Einasto profile for r > 0.01 r_2. 
but for smaller radii deviates to attain a central asymp- 
totic slope —0.86 (see § 4.2 and Fig. 1); in other words, 
a steepening of the density profile is predicted for radii 
r < 0.01 r_2 relative to the flat Einasto shape. In Fig. 11 
we show that a substantially better representation of the 
self-similar solution over the whole range r < r_2 is pro- 
vided by the Sersic-Einasto formula (as recalled in § 1; 
cf. also with Graham et al. 2006) with shape parameters 
T = 0.9 and rj = 0.35. It will be interesting to test such 
a behavior with numerical experiments of higher resolu- 
tions than presently achieved. 

In Fig. 12 we compare the density profile p{r) of the 
self-similar solutions (thick black solid line) to the out- 
comes for six different equilibrium halos extracted from 
the Aquarius A'^— body simulation (Navarro et al. 2010; 
thin colored lines); an Einasto profile with shape param- 
eter 7] = 0.17 is also shown for reference (dashed line). 
The agreement of the self-similar solutions with the sim- 
ulation results is remarkably good for r < a few r_2. A 
discrepancy occurs in the outer regions on approaching 
the caustic discontinuity; this is expected since our fluid 
approximation breaks down there, while the simulated 
halos themselves may be out of eqiiilibrium. 

In Fig. 13 we compare the self-similar profiles (thick 
black solid line) of phase-space density Q{r) to the out- 
comes from the same six Aquarius halos plotted in the 
previous Figure (Ludlow et al. 2010; thin colored lines); 
the powerlaw Q{r) oc r"^'^''^ is also shown for refer- 
ence (dashed line). The self-similar profiles agree with 
the simulation results as to a close powerlaw shape for 
0.04 r_2 <^ a few r-_2- For radii r < 0.01 r_2 that are 
not presently probed by simulations, the self-similar so- 



lutions predict a steepening of Q{r) to attain an asymp- 
totic slope around 2.4 (as expected on the basis of the 
asymptotic analysis of § 4.2, see Fig. 1; cf. also with 
Graham et al. 2006); it will be important to test such 
a behavior with future simulations of higher resolution 
than currently achieved. 

The fact that the powerlaw bcihavior found in simula- 
tions is close to Q{r) cx r^^'*^'' over an extended range 
has often been considered surprising and puzzling; in 
fact, the value —1.875 marks the classic self-similar so- 
lution by Bertschinger (1985) for the purely radial col- 
lapse of a point-mass perturbation. On the other hand, 
the effective spectral index n = 3 corresponding to the 
Bertschinger's solution constitutes a poor approximation 
of the typical overdensity initiating galactic DM halos; 
in addition, the resulting inner density proflle is quite 
steeper than found in galaxy simulations (as recalled in 
§ 1). Here we see that a powerlaw behavior close to 
Q{r) (X r"-'^-®^^ in the halo middle is actually featured 
by our self-similar solutions with more realistic values 
n PS —2.25, when non-radial motions and collision ef- 
fects are included; we find that such a powerlaw behav- 
ior extends over a considerable radial range around r_2- 
Inward of this, Q{r) logarithmically steepens while the 
density asymptotes to slopes < 0.9; at the other end, 
beyond r_2 toward the caustic discontinuity Q{r) fea- 
tures an upturn, for which preliminary evidence has also 
emerged in recent numerical studies (see Ludlow et al. 
2010). 

In Fig. 14 we compare the self-similar anisotropy profile 
(thick black solid line) with the outcomes from the same 
six Aquarius halos plotted in the previous Figures (Lud- 
low et al. 2010; thin colored lines); the outcome obtained 
from the Einasto proflle with rj = 0.17 and the Hansen & 
Moore (2006) /3 — 7 relation is also illustrated for refer- 
ence (dashed line). Once again, the self-similar solutions 
agree with the simulation results (cf. also with the obser- 
vations by Biviano & Poggianti 2009). Note that close to 
the (outer) caustic discontinuity wide oscillations in the 
P = 1 — (Tg/cr^ parameter are found in simulations, as 
expected considering that both the upstream radial and 
tangential velocity dispersions are small. 

5. DISCUSSION AND CONCLUSIONS 

Sharing the widespread drive toward understanding 
the key processes that rule the dark matter (DM) com- 
ponent of cosmic structures from galaxies to their sys- 
tems, we have investigated the detailed structure and 
evolution of cold DM halos. To complement the inten- 
sive and extensive numerical simulations, here we have 
developed advanced models of spherical collapse and ac- 
cretion in an expanding Universe, based on solving time- 
dependent equations for the moments of the phase-space 
distribution function in the fluid approximation. Our 
approach includes non-radial random motions, and most 
importantly, an advanced treatment of dynamical relax- 
ation effects; thus we provide an ciffcctive description of 
the halo dynamical evolution toward equilibrium and of 
the resulting density and velocity structure. Deferring 
to § 5.5 an overall picture of halo structure and develop- 
ment, we discuss first our answers to the hot issues on 
DM dynamics recalled in § 1. 

5.1. Dynamical relaxation 
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Dynamical relaxation involves two mechanisms: phase- 
mixing, due to spreading of neighboring particle orbits 
in phase-space; and violent relaxation, due to irregu- 
lar fluctuations of the gravitational potential. Phase- 
mixing can be traced back to the process of shell cross- 
ing that underlies the classic self-similar treatments (see 
Bertschinger 1985; also Hcnrikscn & Widrow 1999). In 
our fluid approach this is expressed by the left-hand sides 
of Eqs. (11), with action primarily confined to a skin re- 
gion close to the outer caustic; this circumstance enables 
one to describe its effect in terms of boundary jump con- 
ditions of Rankine-Hugoniot type. Violent relaxation, 
instead, is related to collective collisions, expressed by 
the terms on the right-hand side of the third and fourth 
Eqs. (11), and active throughout the volume of the halo. 

The collision terms have been evaluated quantitatively 
with the help of a Fokker-Planck approximation in ve- 
locity space, under conditions of a closely isotropic dis- 
tribution function (see Appendix A2). To derive the 
circularly-smoothed density run p{r), the second Eq. (5) 
shows that it is important to understand the behavior 
of the velocity dispersions cr^ ^i; to that effect, we stress 
that the structure of our Fokker-Planck coefficients en- 
sures erasure of velocity anisotropy after [di {erf. — cr^)]* 
— ((7^ ~ while still conserving the overall random en- 
ergy after (SjCTg)* = — (5t(j,^),/2, see end of § 2.1. 

Anisotropy is initiated at the caustic discontinuity as 
described by the parameter /3+, which specifics the de- 
gree of non-radial DM motions endowed, or acquired 
during the initial infall. The efficiency of their erasure 
is modulated by the strength of the collision terms as 
expressed by the parameter Ki, < 0.1; this reflects the 
amount of gravitational fluctuations induced by dumpi- 
ness in the infalling matter, well beyond the tiny levels 
associated to microscopic graininess. Note that dynam- 
ical relaxation acts like an effective 'torque' mechanism 
after turnaround, in the same vein entertained by Zukin 
& Bertschinger (2010a, 2010b); in fact, these authors flnd 
results qualitatively similar to ours in their parametric 
torque models with decreasing angular momentum and 
vanishing central anisotropy. 

On these grounds, we have found self-similar solutions 
for the spherically-averaged mass density p{r), pseudo 
phase-space density Q{r) = p/a^ and anisotropy param- 
eter /3(r) = 1 — (7^/(7^ • Overall profiles arc obtained on 
matching (as discussed in detail in § 4.4) the self-similar 
solutions corresponding to two different portions of the 
initial DM perturbation, representative of the fast inner 
collapse and of the slow outskirt buildup. We have com- 
pared these overall profiles with the outcomes of state-of- 
the-art A^— body simulations throughout the radial range 
currently probed by the latter, finding a pleasing agree- 
ment. 

5.2. Halo inner regions 

Specifically, in the region inward of a few times r_2 
(the radius where the density slope is —2) we have found 
that efficient dynamical relaxation is the key process to 
produce the following features: a closely universal den- 
sity profile in agreement with the A^— body outcomes, 
well represented in terms of the Einasto formula with 
shape parameter 77 « 0.17; a pseudo phase-space density 
profile with the powerlaw behavior close to Q{r) oc , 



and an anisotropy profile decreasing inward from values 
around 0.25 at r_2 to values /3 « at the center. 

In the very central region r < 0.04 r_2, currently not 
accessible to A^— body simulations, we predict a steepen- 
ing of the mass density profile relative to the flat Einasto 
shape, that would imply a vanishing central slope; the ex- 
pected asymptotic behavior, though with a logarithmi- 
cally slow convergence, reads p{r) oc r""-^^, to the effect 
that the self-similar density profile is best described by 
a Sersic-Einasto formula with shape parameters r = 0.9 
and 77 = 0.35. In parallel, the pseudo phase-space density 
steepens to approach a central behavior Q{r) cx . 
Testing these predictions will require A^— body simula- 
tions of higher resolution than presently achieved. 

5.3. Inner baryonic effects 

Clearly, both our self-similar solutions and A^— body 
simulations refer to pure DM structures. The reader 
must be aware that in the central regions of real galaxies 

small-scale dynamics or energetics related to the astro- 
physics of baryons may alter the baseline DM profiles 
discussed above. These baryonic processes may reconcile 
the outcomes of simulations and self-similar models with 
the observations of galactic dynamics in spiral galaxies 
(especially dwarfs; see Salucci et al. 2007, 2011) and of 
gravitational lensing in galaxies (see Bradac et al. 2009) 
and galaxy systems (see Zitrin et al. 2011; Newman et al. 
2011), which indicate approximately flat density proflles 
already inward of 0.1 r_ 2. 

For example, flattening of the inner density profile may 
be caused by transfer of energy and/or angular momen- 
tum from (baryonic and DM) clumps to the DM field dur- 
ing the galaxy formation process (see El-Zant et al. 2001; 
Ma & Boylan-Kolchin 2004; Tonini et al. 2006). In de- 
tail, upon transfer of tangential random motions from the 
baryons to an initially isotropic DM structure, the origi- 
nal density profile p oc r~^ in the inner region is expected 
to be modified into p cx r'^ [''+2 (2-7) 0o]/[2+(2-7) M . thus 
the final profile is flattened for (3o < 0, down to the point 
of developing a core for ^0 ^ (2 — 7) « —0.3. How- 
ever, a reliable assessment of the amount of angular mo- 
mentum transferred from the baryons to the DM is still 
wanting, and would require aimed numerical simulations. 

Less agreed processes may affect galactic scales r < 10^ 
pc; for example, at the formation of a spheroid, central 
starbursts and accretion onto a supermassive black hole 
may easily discharge enoiigh energy (^ lO®^ erg for a 
black hole mass Af, ~ 10^ Af©) with sufficient coupling 
(> 5%) to blow most of the gaseous baryonic mass within 
r„2 out of the inner gravitational well. This will cause 
an expansion of the DM and of the stellar distributions 
(see Fan et al. 2008, 2010; Ragone-Figueroa & Granato 
2011), so as to flatten the central DM slope. In addi- 
tion, binary black hole dynamics following a substantial 
merger may eject on longer timescales formed stars from 
radii r w 10 (M,/10^ A'/q)^/^^^'^' pc containing an over- 
all mass of a few times the black hole's, and so may cause 
the light deficit observed in some galaxy cores (see Gra- 
ham 2004; Lauer et al. 2007; Kormendy et al. 2009). 

On the other hand, some steepening of the inner den- 
sity profile may be induced by any 'adiabatic' contrac- 
tion of the diffuse star-forming baryons into a disc-like 
structure, as proposed by Blumenthal et al. (1986) and 
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Mo, Mao & White (1998). On the basis of the standard 
treatments, it is easily shown that in the inner region an 
initial powerlaw p oc r~'^ is turned into p oc r"^/*^''"'*'^. 
However, recent numerical experiments (see Abadi et al. 
2010) suggest that the classic treatment of adiabatic con- 
traction is likely to be extreme; actually, in the inner re- 
gion the contraction may be inefficient and the density 
slope hardly modified. 

All these issues are beyond the scope of the present 
paper. 

5.4. Halo outskirts 

At the other end, in the region outward of a few times 
r_2 we find that self-similar solutions with larger e > 1/2 
and less efficient dynamical relaxation apply, consistent 
with the outcomes of numerical simulations. This is ex- 
pected, since the outskirts are built up at later times 
by smoother (i.e., less clumpy) accretion fed on the ta- 
pering wings of a DM perturbation. We note that, as 
highlighted by the scatter in the iV— body results, the 
outskirts structure is subject to a large variance, related 
to the detailed growth histories and to environmental 
conditions wherefrom the infall takes place. 

Relatedly, for redshifts z < 0.5 the accelerating cos- 
mology slows down the time dependence D{t) cx t'^ of 
the growth factor from (i = 2/3to(i«l/2. A strictly 
self-similar solution cannot be obtained in such condi- 
tions, but the overall trend can be captured from not- 
ing that the accretion rate scales as M/M = d/et; thus 
a lower d is equivalent to higher effective values of e, 
which result in steeper density profiles into the out- 
skirts. Moreover, in a flat accelerating Universe, shells 
will be able to turnaround and collapse only if their ini- 
tial position ri is inside a critical radius ta defined by 
5M/M{< ta) = 3/2 (2 Oa)^/^ in terms of the dark en- 
ergy density parameter il^ « 0.7; Subramanian et al. 
(2000) have shown that the resulting outer density pro- 
file is considerably steep, featuring a cutoff toward ta. 

Our self-similar solutions concur with numerical sim- 
ulations in providing a firm basis and in clarifying the 
limitations for the simple dynamical models based on the 
Jeans equation and on a consistent powerlaw shape for 
Q{r) oc r~^'^, see § 1 for details. We find that such mod- 
els are reliable in the significant range 0.01 r_2 ^ & 
few r_2; outside of this, in the outskirts Q{r) flattens ap- 
preciably and is liable to cosmogonical variance, while in 
the inner regions it steepens, though logarithmically. It 
may be interesting to introduce in the dynamical models 
this articulated behavior of Q{r), to investigate how the 
Jeans profiles are affected (taking up from the work by 
Lapi & Cavaliere 2009b). 

5.5. An overall picture of halo structure and 
development 

Finally, we stress the link between the halo structure 
and its two-stage growth history from the vantage point 
of our self-similar solutions, to provide the following over- 
all picture (illustrated in the schematics of Fig. 15). 

The halo formation history starts up with the fast col- 
lapse of a few merging clumps, that simulations indicate 



to be the main contributors to the inner halo mass (see 
Fakhouri et al. 2010; Genel et al. 2010; Wang et al. 
2011); such a strong dumpiness enforces efficient dynam- 
ical relaxation, the key process to produce an approxi- 
mately universal shape of the inner density profile, and 
an isotropic inner halo structure. 

As discussed in § 3, the typical halo mass and radius 
scale as M cx t^/^' and R cx i2/3+2/9e terms of the 
perturbation shape parameter e, or of the corresponding 
effective index of the perturbation spectrum n w 6 e — 3. 
Thus the fast collapse of the inner region affected by 
violent relaxation is seen to correspond to e < 1/6 or 
n < — 2, since for these values the depth of the gravi- 
tational wefl expressed by Vc ~ {GM/Ry/'^ oc f2/9e-i/3 
grows faster than the cosmic time. These values are also 
consistent with the asymptotic analysis of § 4.2. We also 
note that for n < 1 or e < 2/3 the typical specific en- 
ergy v1 cx M'^/^~'^ increases with mass, so that the inner 
region provides an environment conducive to eventual 
melting of the accreted massive clumps. This validates 
spherically-averaged densities, and is consistent with the 
findings from numerical simulations that the residual 
mass fraction still locked into clumps is limited in the 
inner regions (see Springel et al. 2008, their Figs. 11-12). 
Note that the melting of clumps due to tidal disruption 
and dynamical friction may be even more efficient in a 
warm DM framework (e.g., Colin et al. 2000). 

Subsequently, a stage of calmer and smoother accretion 
ensues; this is slower relative to the cosmic expansion 
when n>lore>2/3, since the growth rate scales as 
M/t cx t^/^'^~^. During this stage, simulations indicate 
that most of the mass is gathered about equally from 
minor mergers of many small clumps and from truly dif- 
fuse accretion (see above references). Then the outskirts 
build up from the inside-out, while the preformed inner 
region is nearly unaffected; this is because the smooth- 
ness of the accretion makes dynamical relaxation inef- 
ficient and only mild phase-mixing occurs at the obiter 
caustic (see § 3 and 4.2). Meanwhile, cosmological vari- 
ance and non-radial motions related to the halos' speciflc 
growth history or environmental conditions become rel- 
evant in shaping the halo outskirts, in ways expected to 
be far from universal. 

We envisage such a two-stage framework of halo growth 
to be of much relevance for the astrophysics of galaxy 
clusters, especially in their outskirts as discussed by Lapi 
et al. (2010), and for galaxy/star formation theories; 
e.g., in the latter context it is at the hearth of the two- 
phase galaxy formation picture first proposed by Cook 
et al. (2009). This envisages that the two modes of 
halo growth drive two distinct modes for the evolution of 
baryonic matter, with the development of the spheroidal 
component of galaxies taking place mainly during the 
early fast collapse, and that of the disc component during 
the late slow accretion. 
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APPENDIX 
THE COLLISION TERMS 

In this Appendix we evaluate the cohision terms to be inserted in Eqs. (5) and (11). We take up and adapt to the 
present problem the approach developed by Larson (1969, 1970) in the context of stellar clusters. 

The velocity distribution function 

First of all, we specify the form of the velocity distribution function. We introduce spherical polar coordinate in 
velocity space, and write the velocity components as 



Vr — u = V fl , 
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(Al) 

in terms of the random velocity vector u and the cosine /i of the angle between this vector and the radial direction. 

Wc assume that the normalized distribution /(i/, /i) of random velocities is approximately isotropic, and expand it 
in Legendre polynomials up to the second order 

f{u,^l)c^fo{u)p^{^l) + h{u)p^{^l) + h{y)P2{^x) ; (A2) 

the cocfhcicnts /i and /2 arc assumed to be small corrections to the leading term fo{i') provided for definiteness by a 
MaxwcUian-likc distribution. 
The normalization conditions for / require that 

jd?vf=l , jd^yvf = Q (A3) 

hold, with J d?v = 2n du u"^ J^^ dn being the volume element. These imply that the terms in the expansion 
Eqs. (A2) must depend only on the combinations or u n; thus we base on the expressions 

M-) = J^^^ exp Po(m) = 1 , 

fi{u) = ciufo{u) A(m)=M, (A4) 

where ci and C2 are constants to be determined. 

As discussed in § 2, to close the moment equations we are interested in a distribution with zero skewness, which 
gives Ci = 0. On the other hand, C2 can be determined by considering the second-order moments of the distribution, 
i.e., 

ol = j d^u 1/ V V = o-^ + 2 C2 (tS 

(A5) 

since a^^^ = cr^ + 2 = 3 cr^ and — = 3c2 cr^ hold, we find 

C2=3^^. (A6) 

'''tot 

Fokker-Planck approximation 
Now we turn to the form of the collision terms. In the Fokker-Planck approximation, these are written as 

r.-^ {dt /)* = - E + ^ 2Z ^) (^^^ 

in terms of 

(A8) 

g = 2 J dV fiu')\i.-u'\ . 

The quantity T^, = AwG^mp log A represents the strength parameter of the collective collisions, and depends on the 

clump mass m as m log M/m; here the standard 'Coulomb' logarithm log A has been conveniently written in terms of 
the ratio M/m between the halo and the clump masses (sec Boylan-Kolchin et al. 2008). In fact, numerical simulations 
indicate that the clumps (technically, 'subhalos') in a halo arc distributed close to (see Springel et al. 2008); 

averaging Fj, over such a distribution turns out to be closely equivalent to redefine m as the average clump mass (see 
Ciotti 2010). 
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In axial symmetry, the above terms h and g can be simplified considerably (see Rosenbluth et al. 1956), and to the 
lowest-order one obtains 



3 ~ 47r 



Jo ^ Jv 



(A9) 



correspondingly, the Fokker-Planck equation simplifies to 



Introducing the relevant integrals 



/ 9^9 

V 



+ 



(AlO) 



^0 

Jo 

/>oo 



(All) 



after some manipulation one finds that 



J , 2 



(A12) 



d^g = A'K [X- — + -vK,] ^1.1/5 =—( ^+ 7a ) d^i^i. g = -Sw — . 



J 



J 



Corresponding to the Legendre expansion for / in Eq. (A2), it is convenient to perform the same for (dt /)*, to read 



{dt /). ~ {dt fo). Po(m) + {dt f2h ^2(m) . 

Then the Fokker-Planck equation yields 



(A13) 



r-i(5t/o)*^ 



47r d 
47r d 



/oJ+^(J + z.3/C) 



2,v dv 



(A14) 



u'^ dv 

Taking now the second-order moments leads to the collision terms 



(5ia2), = 87r2r, 



r 

Jo 



dv 



-4(3i/2l+2i/3/C- J)/2 



J 



u{ulC-X)fo-- (51/1-^ ) /; 



-u{uIC-I)fo + - (5ul--] f. 



(A15) 



POO 

{dtal)^ = 8n''r^ / du 
Jo 

we find 



With the adopted form of / discussed in § Al, we find 



(A16) 



so that the collision terms finally read 



^2 „2 



{dt <)^ = -— G^mp log A 3 ^ , 

crtr>t 
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(A17) 



{dt<jl\ = + G^mplogA " 3 ^ 

15 CTiLt 



With the effective number of clumps A/' = M/m approximately constant, the above collision terms do not break 
self-similarity, and can be written as 

^tot 

(A18) 

i^K = +^ ^ in - n) ■ 



in terms of the self-similar variables defined in Eqs. (10) of the main text. The strength parameter of dynamical 
relaxation 



= 64V3^1ogA^ logA^ 
405 A^ A^ 

depends mainly on the dumpiness \/N of the infalling matter. In Fig. Al we illustrate as a function of A/", and 
highlight the regimes corresponding to the early fast collapse and to the late slow accretion. In more detail, during 
the early fast collapse of the halo inner region, a limited effective number of clumps A/^ < 10 is expected to enforce 
efficient dynamical relaxation with k^, « 0.1, while during the late slow accretion dominating the growth of the halo 
outskirts, many small infalling clumps with A/" > 30 make dynamical relaxation much less efficient with < 0.01. 
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Figure 1. Asymptotic inner slopes of density and (pseudo) phase-space density as a function of e, for different values of tlie central 
anisotropy parameter. We recall from § 2.1 that efficient dynamical relaxation enforces /3o = 0. 
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Figure 2. The self-similar solution for e = 1 and /3+ = 1 corresponding to purely radial orbits, in the absence of collisions (re* = 0); this is 
just the case considered by Bertschinger (1985). Green dotted lines ilUustrate the expected asymptotic behaviors (see § 4.2). Vertical red 
lines mark the radius -R200 where the density amounts to 200 times the background's. Same linestyles are adopted in the following figures. 
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Figure 3. The self-similar solution for e = 1/6 and /3+ = 1 corresponding to purely radial orbits, in the absence of collisions (ft* = 0). 
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Figure 4. The self-similar solution for e = 1/6 and /3+ = 0.5 in the absence of collisions re* = 0. In the middle-right panel, solid line 
refers to Q{r) defined in terms of the total velocity dispersion, while dashed line refers to Q{r) defined in terms of the radial dispersion. 
In the bottom-left panel, solid line refers to the total velocity dispersion, dot-dashed line to the radial component and dashed line to the 
tangential component. 



20 



Lapi & Cavaliere 





1.000 



Velocity Dispersions 




0.001 



0.010 0.100 
X 



1.000 



Binney porometer 




0.001 0.010 0.100 1.000 



Figure 5. Self-similar solution for e = 1/6 and fi-^. = 0.25, in the absence of collisions = 0. 
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Figure 6. Self-similar solution for e = 1/6 and fi-^. = 0.25, in presence of collisions with strength parameter = 0.01. 
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Figure 7. Self-similar solution for e = 1/6 and fi-^. = 0.25, in presence of collisions with strength parameter = 0.05. 
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Figure 8. Self-similar solution for e = 1/6 and fi-^. = 0.25, in presence of collisions with strength k* = 0.1. 
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Figure 9. Self-similar solution for e = 1/8 and = 0.25, in presence of collisions with strength k* = 0.1. 
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Figure 10. Position of the caustic discontinuity as a function of e, for two values of l3+ and k*; these parameters only mildly affect the 
caustic position at given e. The curves have been computed at discrete values of e, as highlighted with the open symbols; the specific 
values e = 1/8 and 1/2 used in the next figures are marked by filled dots. Note that for values e > 1/2 the caustic position are close to the 
standard virial radius iJ200 ~ fita/2 recalled in § 3. 
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Figure 11. Comparison of the inner self-similar density profile (top) and density slope (bottom) with standard fitting functions to the 
equilibrium outcomes of numerical simulations. Solid line refers to the self-similar solution with parameters e = 1/8, /9+ = 0.25, = 0.1, 
dot-dashed line refers to the NFW profile, dashed line refers to the Einasto profile with shape parameter -q = 0.17, dotted line refers to the 
Sersic-Einasto profile with shape parameters r = 0.9 and r) = 0.35. The shaded area highlights the radial range not accessible to present 
numerical simulations. 
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Figure 12. Comparison of the self-similar density profiles (thick solid line) with e = 1/8, (5+ = 0.25, ft* = 0.1 (region inward of r—2) and 
with e = 1/2, /3+ = 0.25, k* = 0.01 (region outward of r—2) to the outcomes for six different halos extracted from the Aquarius Af— body 
simulation (Navarro et al. 2010; thin colored lines); the standard Einasto profile with rj = 0.17 is also illustrated for reference (dashed line). 



28 



Lapi & Cavaliere 

Pseudo phase-space density 



4 



O 



a 

O 







2 - 



■4 



^ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 

- 


This work 




- 

- 
- 


_ Powerlaw 1.875 : 




s 


N-body (Ludlow+ 11) 




s 

s, 


1 1 1 





N _ 



0.01 



0.10 1.00 
r/r_2 



10.00 



Figure 13. Comparison of the self-similar profiles of pseudo phase-space density (thick solid line) with e = 1/8, /9+ = 0.25, = 0.1 
(region inward of r_2) and with e = 1/2, /9+ = 0.25, k* = 0.01 (region outward of r_2) to the outcomes for six different halos extracted 
from the Aquarius TV— body simulation (Navarro et al. 2010; thin colored lines); the powerlaw Q oc with x = 1.875 is also illustrated 
for reference (dashed line). 
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Figure 14. Comparison of ttie self-similar anisotropy profiles (thick solid line) with e = 1/8, /?+ = 0.25, k* = 0.1 (region inward of r_2) 
and with e = 1/2, /9+ = 0.25, re* = 0.01 (region outward of r_2) to the outcomes for six different halos extracted from the Aquarius 
A'^—body simulation (Navarro et al. 2010; thin colored lines); the anisotropy profile obtained on combining the Einasto density profile with 
7y = 0.17 and the Hansen & Moore (2006) /? — 7 relation is also illustrated for reference (dashed line). 
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Figure 15. Schematics of the link between the halo structure and its two-stage growth history as it emerges from our self-similar solutions 
and from numerical simulations (see § 5 for details); note that cosmic time runs from left to right. 




Figure Al. The strength parameter of the dynamical relaxation as a function of the effective number of clumps M = 
matter, see Appendix A for details. 
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